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ABSTRACT 

Almost all hydrodynamic accretion disk models parametrize viscosity with the dimen- 
sionless parameter a. There is no detailed model for a, so it is usually taken to be a constant. 
However, global simulations of magnetohydrodynamic disks find that a varies with distance 
from the central object. Also, Newtonian simulations tend to find smaller a's than general 
relativistic simulations. We seek a one-dimensional model for a that can reproduce these two 
observations. We are guided by data from six general relativistic magnetohydrodynamic ac- 
cretion disk simulations. The variation of a in the inner, laminar regions of the flow results 
from stretching of mean magnetic field lines by the flow. The variation of a in the outer, 
turbulent regions results from the dependence of the magnetorotational instability on the di- 
mensionless shear rate. We give a one-dimensional prescription for a{r) that captures these 
two effects and reproduces the radial variation of a observed in the simulations. For thin disks, 
the prescription simplifies to the formula a{r) = 0.025 [q(r)/ 1 .5] , where the shear parameter, 
q(r), is an analytical function of radius in the Kerr metric. The coefficient and exponent are 
inferred from our simulations and will change as better simulation data becomes available. 
We conclude that the a-viscosity prescription can be extended to the radially varying a's ob- 
served in simulations. It is possible that Newtonian simulations find smaller ar's than general 
relativistic simulations because the shear parameter is lower in Newtonian flows. 

Key words: accretion, accretion discs, black hole physics, hydrodynamics, (magnetohydro- 
dynamics) MHD, gravitation 



£N| . 1 INTRODUCTION 



Accretion disk magnetic fields and turbulence act as a large scale 
viscosity, draining angular momentum and energy from accreting 
gas. The magnitude of this viscosity is uncertain, w hich adds a free 
param eter to accretion disk models. For example, IPringle & Reed 
(1972) leave the ratio of the disk's radial velocity to its circular 
velocity as a free parameter. They call this ratio y/100 and esti- 
mate y ~ 1. They note that y might be radius dependent and that 
a strong radial dependence could change qualitative features of the 
disk. However, other important features of the disk, such as its lu- 
minosity, depend weakly or not at all on y, so progress is possible 
with out a detailed model fory. 

IShakura & Sunvaevl jl973h parametrized the viscosity with a, 
the ratio of stress to pressure. It is related to Pringle and Rees's y 
by y = I00a(h/r) 2 , where h/r is the disk opening angle. 0They an- 
ticipated that a will be a function of radius, citing experiments of 
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turbulence i n Taylor-Coutte flow, flow between rotating cylinders 
jTavloi|[l936h . These experiments show that the torque exerted by 
a turbulent flow on the cylinders depends on the rotation rate of the 
cylinders and the separation between them. The torque in the exper- 
iment is related to a, so a should depend on accretion flow prop- 
erties that vary with radius. The analogy is too rough to produce a 
quantitative model, so they assume a is a constant for simplicity. 

The most significant theoretical breakthrough towards an un- 
derstanding of a has been the realization that the magnetoro- 
tational instabilit y (MRI) can drive turbulence in ionized accre- 
tion disks jyelikhovll 19591 : Ichandrasekh"a3ll96d : iBalbus & Hawlevl 
Il99llll99ii) . A weak seed magnetic field and a radially decreasing 
angular velocity profile are all that are required to trigger the insta- 
bility and both are present in disks. This makes it possible to run 
magnetohydrodynamic (MHD) disk simulations without using the 
a-viscosity prescription. 

Nonetheless, because of its relative simplicity, the a- 
viscosity prescription has not dimmed in importance. Over the 
past year, more than 300 papers cited the pioneering work of 
Shakura & Sunyae\] J 1973h and more than 1000 papers made ref- 
erence to some disk solution based on the «-viscosity prescript ion 
(such as the relativistic thin disk of iNovikov & Thornd H973I) or 
the advection dominated accretion flow (ADAF) of Naravan & Yil 
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(1994)). However, despite its interest, there is still no widely ac- 
cepted model for the size or shape of a. In applications, it is 
typically assum ed to be a constant between 0.01 and 0.1 (e.g., 
lGouetal.ll201ll) . 

There are two varieties of MHD simulations, local and 
global, and both provide hints about the size and shape of a. 
A standard setup for a local simulation is a box of weakly 
magnetized fluid in the shearing sheet approx i mation (e.g., 
Hawlev & Balbu sl Il99ll: iBrandenburg et all 1 1 9951; lHawlev et al 



199511 1994 IStone & Balbuslll996l : lBrandenburdl200ll ; ISano et al 



2004). The dimensions of the box are typically a few disk scale 



heights. The earliest global simulations of MRI turbulent disks 
used a Newtonian potential jArmitagd 1 19981 : iMatsumotcl 1 19991 : 
lHawlevl2qodlMa chida et al. 2000) and, later, a pseudo-Newtonian 
potential faawlev & Krolikl 1200 ll) . The first general relativisti c 
global simulations were carried out by |Pe Villiers et alj j2003h . 
A standard setup for global simulations is a torus of fluid in hy- 
drostatic equilibrium, threaded with a weak magnetic field (e.g., 
|Pe Villiers et alj2003l) . In local and global simulations, differential 
rotation triggers the MRI and drives turbulence. When the simula- 
tion reaches a quasi-steady state, one can compute the ratio of stress 
to pressure and so measure a. Local simulations focus on resolv- 
ing small scale physics, such as saturation of the MRI, while global 

simulations attempt a complete portrait of the accretio n disk. 

Global simul ation s have h inted at the shape of a. lPenna et al.1 

j20ld) and lPenna et al. (2012b) presented general relativistic MHD 
(GRMHD) simulations of MRI turbulent disks in the Kerr met- 
ric, the spacetime of a spinning black hole. They c o nfirme d that 
the relativistic cr-disk model of iNovikov & Thornd Jl973h gives 
a good description of the simulation data. This was partly mo- 
tivated by earlier suggestions that magnetic stresses at the inner 
edges of MHD disks invalidate the assum ptions of cr-disk mod- 
els (e.g.. lKrolikli999l : lAgol & Krolikl 200d) . The GRMHD simula- 
tions produced lumin osity inside the innermost stable circular or- 
bit (ISCO), where the INovikov & Thornd d 19731) model is entirely 
dark, but it was a modest contribution to the total luminosit y of the 
disk (on the order o f a few percent; iKulkarni et alj|201 1 ). Later , 
IPenna et al.1 J2012bl) generalized the INovikov & Thornd Jl973l) 
model to self-consistently incorporate nonzero luminosity inside 
the I SCO, bringing it clo ser to the GRMHD simulations. 

IPenna et alj d2010h and IPenna et aD d2012bl) measured the 
shape of a. For accretion onto a non-spinning black hole, they 
found a is small at the event horizon, increases to a maximum near 
the photon orbit, declines to ~ 0.06 near the ISCO, and then contin- 
ues to decline, albeit more slowly (see, e.g., Figure 7 of IPenna et"al] 
l2012bl) . In a com pletely different contex t, MHD simulations of pro- 
toplanetary disks. iFromang et all d20 1 lb found a radially varying a 
with similar shape. However, the overall size of their a was over an 
order of magnitude lower, peaking at 0.013 and declining to below 
0.002. 

The val ue of a is notoriously difficult to pin down 

( IPessah et all l2007h . Local simulations find that a is strongly 
affected by grid sca le di ssip ation as well as stra tifi cation 
l Lesur & LongarettH |2007 | : Fromang & Papaloizoul 120071 : 
ISimon & Hawlevl 120091 ; iDavis et al.1 |201C|) . In global and lo- 



cal simulations, a depends non-monotonically o n resolution at all 
but the highest resolutions dSorathia et alj|2012l) . In models with 
net magnetic flux, a scales with increas ing flux ( lHawlev et al.1 
1 19951 ; ISano et alj|2004l ; IPessah et alj|2007h . Finally, or depends on 
the in itial magnetic field geometry and strength ("e.g. JSorathia et all 

A combination of these effects can probably explain some 



of the discr epancy between the a values found by IPenna et"al] 
d2012bl) and Fromang et al. U2011I) . The former employed an ini- 
tially poloidal magnetic field with initial gas-to-magnetic pressure 
ratio p = 100 and resolution 256 x 64 x 32 in (r,0,0). The later 
employed an initially toroidal field with initial /? = 25 and resolu- 
tion 512 x 256 x 256. Neither had explicit dissipation or net mag- 
netic flux. Initially toroidal fields tend to beget smaller o-'s than 
initially poloidal fields, but a also tends to be inversely related to 
the initial p. The sc aling of a with resolution is non-monotonic 
jSorathiaet"aT]|2012h . It is not clear whether these effects are large 
enough to explain why the a's measured from the general relativis- 
tic simulations are over an order of magnitude larger than the a's 
measured from the Newtonian simulations. 

In this paper, we present a one-dimensional model for the 
shape of a and we show that relativistic corrections enhance the 
a's measured in GRMHD simulations relative to Newtonian simu- 
lations. The one-dimensional model has two components. The first 
component is generated by large-scale, mean magnetic fields and 
is based on the model of lGammid ( 119991) . It dominates in the inner 
regions of accretion flows, where plunging gas stretches and am- 
plifies the frozen-in magnetic field. The second, "turbulent," com- 
ponent describes the dependence of a on the shear rate of the flow. 
The dependence of a on the shea r rate in turbulent flows wa s ear- 
lier emphasized bylGodonl d 1995T ). lAbramowicz et alj ( 119961) , and 
IPessah et all d2008l) . Combining the mean field and turbulent com- 
ponents yields a one-dimensional model for the shape of a. 

We use this model to fit the profiles of a versus radius ex- 
tracted from six GRMHD simulations. The simulations describe 
thin disks accreting onto a non-spinning black hole at two differ- 
ent resolutions, a thin disk and a thick disk accreting onto a spin- 
ning black hole (with dimensionless spin parameter a/M = 0.7), 
and thick disks with two different initial magnetic field topologies 
accreting onto a nonspinning black hole. We show how the radial 
variation of a changes the structure of tr-disk solutions. Finally, we 
note that the enhancement of a seen in GRMHD simulations rela- 
tive to Newtonian simulations can be explained by the dependence 
of a on shear rate. 

There is another reason a is interesting which we have not 
yet mentioned. Turbulent fluids can be complicated. They are dis- 
ordered solutions of nonlinear equations that require great effort to 
solve numerically. And yet, from fully developed turbulence, sim- 
ple scaling laws emerge with apparently universal properties. For 
example, shearing box simulations of MRI-driven turbulence with 
a Keplerian rotation profile always find that the ratio of Maxwe ll 
stress to Reynolds stress is a constant, « 4 ( IPessah et al.l l2006b). 
The same ratio appears independently of the magnitude or geom- 
etry of the magnetic field. It see ms only to depend (in a simple 
way) o n the shear rate of the flow ( lHawlev et aljl 999; Pessa h et al] 
2006b). Similarly, the viscosity parameter, a, obeys a remarkable 
scalin g law: a(l ~ 1/2, where is the gas-to - magnetic pressure ra - 
tio jBlackman et al]|2008l ; lGuan et al J l2009l ; ISorathia et alj|2012l) . 
Formulas this simple should have simple explanations; they should 
not just appear at the end of large numerical calculations, as if by 
coincidence. We hope that clarifying some of the physics underly- 
ing a will improve our understanding of these scaling laws. 

The paper is organized as follows. In ^2] we give an overview 
of physics in the Kerr metric. In i}3] we describe our one- 
dimensional model for a(j). The GRMHD simulations are de- 
scribed in § 34161 we give a broad overview of the simulations in 
^4] analyze the two fiducial simulations in ^5] and analyze the re- 
maining four simulations in jj6] In Sj7] we show how a radially vary- 
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ing a affects the structure of a-disk solutions. We conclude with a 
summary and discussion in ^8] 



2 PRELIMINARIES 

For our investigation of a, we will need to compute the angular ve- 
locity, shear rate, and epicyclic frequency of an accreting gas in the 
Kerr metric, so we review their definition here. This also helps to 
establish notation. As an intermediate step, we discuss the transfor- 
mation between the Boyer-Lindquist and fluid frames. 
The Kerr metric in Boyer-Lindquist coordinates is 

ds 2 = - (1 - IMrjT) dt 2 - (4Mar sin 2 dtd(f> + (2/A) dr 2 
+ (r 2 + a 2 + 2Ma 2 r sin 2 sin 2 Od<jr . (1) 

Here M is the mass of the black hole, a is its angular momentum per 
unit mass (0 < a < M), and the functions A, X, and A are defined 
by 



where, 



A ee r 2 - 2Mr + a 2 , 

2 ee r + a cos - 9, 

A ee (r 2 + a 2 f - a 2 A sin 2 6. 



(2) 
(3) 
(4) 



The accreting, magnetized gas is characterized by its four- velocity, 
ur 1 , density, p, pressure, p, and internal energy, u, and by the elec- 
tromagnetic field, F 1 ™ . We set G = c = 1. 
The angular velocity of the gas is 



dt u> 



(5) 



Circular equatorial geodesies have CI = M [/2 /(r 3/2 + aM' /2 ), and 
the motion of a geometrically thin disk is well approximated by 
circular geodesies outside the ISCO. However, motion inside the 
ISCO and the motion of thick disks are not Keplerian, so we leave 
CI unspecified for now. 



2.1 Inertial Fluid Frame 

We would like to evaluate the shear rate and epicyclic frequency in 
the inertial fluid frame, rather than the Boyer-Lindquist, ZAMO, or 
any other frame, because the physics is simplest there. Also, this is 
the frame where a is defined. In this frame, the equivalence prin- 
ciple lets us ignore gravitational forces at the center of the fluid 
parcel we are following. Tidal forces and other gravitational forces 
that become important over large distances will not concern us be- 
cause shear rate and epicyclic frequency are local measurements. 

Measurements in the Boyer-Lindquist frame, (dt, dr, dS, dip), 
are related to measurements in the inertial fluid frame, 
{oJ, of, of, o>^\ by the the transformation matrix oJi jKrolik et al.l 
lioOSt iBeckwith et aU 2008 ; Kul karni et alj201lh : 



of. = Im ,u ,u , w\ , 

oJi = — iyUyU 1 , 1 + u' U r + U Ug, 0, U r U?} , 
(Ji = — (ugll' , UgU' , 1 + UgU 9 , Ugll^ , 



-5 = ^- 

* N 3 



(-{,0,0,1), 



(6) 
(7) 

(8) 

(9) 



N 



s = -Co/ \C \ , 

grr^gltC 2 +g rr Cl + 

N 3 = ^gnt 2 - 2g t<1 { - 



I = Uf/U,, 



+ 2gyCiC z , C = u'u, + u^u^, 
Ci = u'u,, 
C2 = u r U{,. (10) 



Hatted indices refer to fluid frame quantities and unhatted indices 
refer to Boyer-Lindquist frame quantities. In the orthonormal fluid 
frame the metric is the Minkowski metric, 7/ a g = diag(-l, 1,1,1). 
So hatted indices are raised and lowered with the Minkowski metric 
and unhatted indices are raised and lowered with the Kerr metric. 

The fluid frame basis ©-(O was constructed using a Gram- 
Schmidt process. There is some arbitrariness in the orientation of 
the frame, but we followed standard conventions. Equation ([SJ for 
oj; is necessary because the Lorentz factor in the fluid frame should 
satisfy 7 = -;(,- = -cJiu^ = 1. The next step in the Gram-Schmidt 
process is to define oj^ such that it is orthogonal to ojj and has 
no component along dr or dO. Finally, oj f and a> g are constructed. 
Notice that oj- r has a nonzero component along d(f>, and oj e has 
components along all four Boyer-Lindquist directions. This is un- 
avoidable. However, the most important directions for the physics, 



and 



are aligned as closely as possible with their Boyer- 



Lindquist analogues. 

The arbitrariness in the construction of the fluid frame leads to 
an ambiguity in the definition of a. One usually avoids quantities 
with these sorts of ambiguities. But, as discussed in ^T] a is too 
useful for accretion disk modeling and turbulence theory to aban- 
don. So our strategy is to define a in the most natural way possible 
and see where this leads. 

If the poloidal velocity is much smaller than the az- 
imuthal velocity, as usually happens everywhere except in the 
innermost regions of the disk, the fluid frame basis simplifies 
dNovikov & Thornelll973l) : 



OX = u,dt + u,pd(p, 
0/ = D- 1/2 dr, 
J = rd6, 

J = yrJl 1 ' 2 (d(/> - Qdt) . 



(11) 
(12) 
(13) 
(14) 



The Lorentz factor is 7 = -\J-g,,u' and 

31 = 1 +aV 2 + 2MaV 3 , D = 1 -2Mr l +aV 2 . (15) 

The relativistic factors 3K, T), and 7 are unity at large radii. In fact, 
in this limit, the fluid frame basis is exactly aligned with the Boyer- 
Lindquist frame. So the ambiguities in the general relativistic defi- 
nition of a discussed earlier are not important in the outer regions 
of the disk. There remains a potential ambiguity in the construction 
of the Boyer-Lindquist r and 6 coordinates, even when the gas is 
nonrelativistic, because the black hole's spin axis breaks the spher- 
ical symmetry of spacetime. But this is also negligible far from the 
black hole where frame dragging is weak. 



2.2 Electric and Magnetic Fields 

All observers interact with the same electromagnetic field. They 
might measure its components, F^ v , differently, depending on their 
reference frames, but the underlying object is the same multilinear 
map, F. 

The electric and magnetic fields are not so universal: each 
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observer splits the electromagnetic field into different electric and 
magnetic components. An observer with four-velocity if measures 
electric and magnetic fields 



UyF"' 



If = u v ,F"\ 



(16) 



where the dual Faraday tensor is ,F' ly = € t ""' A F KA /2, the Levi-Civita 
symbol is e*"'*' 1 = - [/jwc/1] / s^g, and [juvk/1] is the completely 
antisymmetric symbol, equal to either 0, -1, or +1. 

The electric and magnetic fields ef and V transform as ten- 
sors, so they can be evaluated in any frame. In general, all four 
components are nonzero. The four- velocity of the observer, if, 
appears in J16t . so each observer interacts with different electric 
and magnetic fields. Following standard convention, we denote the 
Boyer-Lindquist observer's splitting W and W, and we denote the 
fluid frame's splitting and If. In Boyer-Lindquist coordinates 
B' = E = 0, and in the fluid frame e' = b' = 0, by the antisymme- 
try of F MV . The remaining nonzero components in these frames are 
the usual three- vectors of special relativity. 

In ideal MHD, e p = 0, so the fluid frame splitting of F^. is usu- 
ally the simplest. For calculations, Boyer-Lindquist coordinates are 
sometimes more convenient than the fluid frame. So it is common 
to work with If, the fluid frame observer's magnetic field expressed 
in Boyer-Lindquist coordinates. This is legitimate, because If is a 
4-vector, but no observer would measure b' , V , b e , or b*. Physically 
meaningful quantities are W ,b e , and If. 

It is possible to convert between B M and b" directly, without 
reference to F 1 ™. In Boyer-Lindquist coordinates, the formulae are: 



b' = BVg, v , V 
B" = Wu' - b'if. 



B> + b'li 



(17) 
(18) 



In these equations, i runs over r, 8, <f>, and ji, v run over t, r, 8, (f>. 



frame, u^(r) = 0, so the derivative is0: 



1 u*(r + df) 
- lim . 

2 m dr 



(21) 



We use equations dl lt-dl4b to rewrite fluid frame measurements in 
terms of Boyer-Lindquist measurements, obtaining: 



1 



-y 2 WrQ. r . 



(22) 



This is the product of th e Newtonian shear r a te rfl , r /2 with a rel- 
ativistic correction, y 2 JK. iNovikov & Thome] i 19731) state this for- 
mula without proof. It describes the rate of change of a disk's angu- 
lar velocity with radius, assuming the azimuthal velocity dominates 
the poloidal velocity. 

A dimensionless measure of the shear rate is 



-2oW£2 : 



, d log fi 



(23) 



d log r 

Positive q corresponds to angular velocity decreasing with radius. 
Solid body rotation is q = 0. Accretion disks generally have posi- 
tive q, although q changes sign near the photon orbit in black hole 
disks . Flows with positive q are unstable to the M RI fVelikhov 
ll959l : IChandrasekhaj|l96(]|:lBalbus & Hawlevlll99ll) , and the value 
of a is a function of a jPessah et alj|2008l) . Flows with q > 2 
are Rayleigh (hydrodynamic ally) un s table. The MRI has been ana- 
lyzed in the q > 2 regime bv lBalbuslJifJl 2]j. Circular geodesies in 
Newtonian gravity ha ve q = 3/2. Ci rcular, equatorial geodesies in 
the Kerr metric have dGammiell2004l) 



1 = o 



2 1 - 3Mr~ l + 2aM l ' 2 r-3l 2 ' 



(24) 



which becomes 3/2 at large radii. The fact that a is a function of q, 
and in the Kerr metric q is a function of radius, is the basis for the 
turbulent component of the one-dimensional a model in ^3] 



2.3 Shear Rate 

The shear rate is a measure of how rapidly the angular velocity 
of an accretion disk varies with radius. It is rQ. r /2 in Newtonian 
gravity. More generally, one can define the special relativistic shear 
tensor 




where h a p = g a p + u a Us is the projection ten sor and © = u a - a is 
the expansion scalar ( Nov ikov & Thornell973r) . Then the shear rate 
is the r0 component of this tensor measured in the fluid frame: 
(Tfi. The shear tensor so defined is trace-free and symmetric. To 
obtain the general relativistic version, one should replace partial 
derivatives in equation ( 119t with covariant derivatives. We do not 
require this generality because in the inertial fluid frame the laws 
of physics take on their special relativistic form without gravity, by 
the equivalence principle. 

To obtain a simple formula for the shear rate, let us assume 
the poloidal velocity is small and the flow is axisymmetric. Then 
we can ignore the expansion scalar, 0, use equations il lt-dl4t for 
the fluid frame transformation, and ignore derivatives with respect 
to The fluid frame projection tensor is hgjj = diag (0, 1, 1, 1), so 
the shear rate is 

<r* = (20) 
which resembles the Newtonian shear rate, rSl r /2. In the fluid 



2.4 Radial Epicyclic Frequency 

The radial epicyclic frequency is the frequency at which a radi- 
ally displaced fluid pa rcel will oscillate. As a function of q, it is 
jPringle & Krn3l2007h 

k= V2(2 - q)£i. (25) 

This expression is valid in the inertial fluid frame provided we use 
the relativistic expressions for CI and q, equations lO and {23}. 

Circular geodesies at the ISCO have k = by definition, which 
makes it easy to see that they have q = 2.A thin accretion disk starts 
with q = 3/2 at large radii and then q increases until it reaches 2 at 
the inner edge. T he radial depend ence of k 2 for circular, equatorial 
Kerr geodesies is (Gammie 2004) 

2 _ 1 l-6M/r + 8aM"V 3/2 -3a 2 /r 2 
K ~r> 1 - 3Mr' + 2aM l l 2 r - 3 ' 2 ' ( ' 

The epicyclic frequency is zero at the ISCO, and imaginary inside 
the ISCO, signaling the instability of circular geodesies. In ^5]and 
^6] we compute the epicyclic frequency of simulated GRMHD ac- 
cretion disks as a function of radius. The inner edge of the disk 

2 We abuse notation for clarity. The commutators of fluid frame basis ele- 
ments do not vanish in general, so they are not coordinate induced: there is 
no "r" coordinate satisfying df = e , 

3 Note that our simulated disks have q < 2 at all radii. The scaled epicyclic 
frequency reaches a minimum near the ISCO and increases in the plunging 
region. 
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can be identified with the minimum k. Only for thin disks does k 
have a sharp minimum at the ISCO. In general, pressure gradient 
forces and magnetic stresses smear out the inner edge of the disk 
and displace it from the ISCO. 



3 ONE-DIMENSIONAL MODEL FOR a(r) 

In this section we define a one-dimensional prescription for the de- 
pendence of a on radius. Our model is the sum of two components: 
a turbulent component that dominates in the outer regions of the 
disk and a large scale magnetic field component that dominates in 
the inner regions of the disk. We discuss each component separately 
before combining them into a single prescription for a(r). We will 
fit this a(r) prescription to data from GRMHD simulations in §ij4} 



6 



The standard a viscosity prescription is 



: ap, 



(27) 



where Tjx is the fluid frame stre ss, p is pressure, and a is a con- 
stant (Shak ura & Sunvaev|[T973h . There are two equivalent ways 
to modify this prescription. One can add extra factors multiply- 
ing the RHS of equation ( 1271 ) and keep a a constant, or one can 
keep equation ((27} unchanged but define a to be a function of ra- 
dius, a = a(r). We chose the later convention because it makes it 
easy to adapt old a-disk solutions to the new prescription: just in- 
sert a(r) wherever a appears. Both approaches have been used in 
the past so it requires a bit o f care to compare results. For exam- 
ple, [Abramowicz et al. ( 1996) discussed the dependence of a on q, 
while Pes sah et"al] {2008) modified j27b and kept a a constant. 



3.1 Turbulent a 

We have seen in 32.31 that q can be a function of radius, either 
through relativistic corrections to the Newtonian shear rate, or 
through the dependence of CI on radius. It turns out a is a func- 
tion of q, and so it too can depend on radius. That a depends on q 
is perhaps not surprising. MHD flows are unstable to the MRI and 
have finite a when q > 0, but they are stable against the MRI and 
have vanishing a when q < 0. Some dependence of a on q must 
conn ect these two regim es. 

IPessah et al.1 {2008) examined the dependence of a on q nu- 
merically. They used nonrelativistic shearing box simulations with 
resolution 32 x 192 x 32 in r x <p x z, zero net magnetic flux, and 
initial gas-to-magnetic pressure ratio 8 = 200. They ran a series of 
simulations in which they varied q from q = —1.9 up to q = 1.9 in 
steps of Aq = 0.1. Each simulation ran 150 orbits and they mea- 
sured a by averaging the data from the last 100 orbits. For q > 0, 
they found the power-law scaling 



a cc q 



(28) 



where n is between 2 and 8. Higher resolution simulations are 
needed to determine n more precisely. For q < 0, they found a = 0, 
as expected from the MRI stability criterion. 

Between the inner edge of a thin accretion disk and the outer, 
non-relativistic regions, q only varies by 50% (c.f. 92.3b . but a can 
vary by much more if the exponent in equation l |28t is large. For 
n = 8, the change in or is a factor of 10. 

To get a quantitative prescription for a(r), we need q(r). So 
first we solve standard a-disk equations with constant a = a . This 
gives q(r). Then we define a(r) = or [q(r)/1.5] n . The exponent is 
a free parameter to be determined from MHD simulations. Now 



one could iterate: feed a(r) back into the a-disk equations, and re- 
evaluate q(r). For simplicity, we do not iterate. We check our model 
for a(r) against GRMHD simulation data in § 94161 but the simula- 
tion data are too noisy to justify computing a(r) more precisely for 
now. 

The a-disk solutions we use ar e the relativ i stic sl im disk so- 
lutions of Abramowic z et al. | dl988l) : ISadowskil d201ll) . This is a 
family of one-dimensional solutions for black hole accretion with 
four free parameters: black hole mass and spin, accretion rate, and 
a. At low accretion rates t hey r educe to the standard thin disks of 
IShakura & Sunvaevl d 19731) and lNovikov & Thornd dl973l) . and at 
high accretion rat es they become advection dominated and are sim- 
ilar to slim disks (Abra mowicz etaljl9 88: Sadowski 2oTl|) 

This prescription for a(r) neglects the contribution of large 
scale, mean magnetic fields, which can exist even in laminar flow. 
These are important near and inside the inner edge of the disk, 
where the flow acquires a large radial velocity and stretches the 
frozen-in magnetic field. We discuss this contribution to a(r) in the 
next section. 



3.2 Mean Magnetic Field Stresses 

Penn a et al. I d20ld) observed large scale, mean field stresses in 
GRMHD simulations of black hole accretion disks, and they 
showe d that a one dimensional model developed by iGammiel 
dl999h could fit the d ata. 

iGammid d 19991) solved for the motion of a fluid with a frozen- 
in magnetic field as it plunges into a Kerr black hole along the equa- 
torial plane. The inner boundary of the flow is at the event horizon 
and the outer boundary of the flow is at r B , where the flow is as- 
sumed to have zero radial velocity and Keplerian angular velocity. 
The governing equations are mass, angular momentum, and energy 
conservation, and Maxwell's equations. There is no dissipation and 
the pressure and internal energy of the gas are neglected. The solu- 
tions are time-independent, axisymmetric, and vertically averaged. 
They provide the rest mass density, velocities, and magnetic field 
of the flow as a function of radius. The free parameters are a/M, r B , 
and the amount of magnetic flux threading the horizon. Following 
IPenna et al.l d20ld) . we parametrize the flux threading the horizon 
by 



/ \B r \dA 



(29) 



which is dimensionless. A "magnetically arrested disk" cor- 
responds to 0bh = V^T >. 50 dNaravan et al] 120031 : 
iTchekhovskov et al.ll201 ll:lNaravan et alj2012h . 

Gammie's code for generating solutions numerically using the 
shooting method is available on the web. The solutions are ex- 
pressed in terms of B', but f follows from equations l|17|l-l|18t, and 
the stress, —bfbx, follows from equations l[6]l-l[9}. 



3.3 Combined Model 

Combining the turbulent and mean field contributions to a(r) gives 
the one-dimensional prescription 



a(r) = a 



q(r) 



3/2 



- ai- 



b f (r)b t (r) 
p(rY ' 



(q > 0). 



(30) 



http: //rainman. astro, illinois . edu/codelib/codes/inflow/src/ 
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If there are no large scale fields and q = 3/2, then a(r) = a , a 
constant. If q < 0, then one sho uld only include the second term 
on the RHS. The lGammid dl999l) solutions do not include gas pres- 
sure so we have divided the mean field stress by p(r) r , which is 
proportional to pressure for a polytropic gas. This is a crude substi- 
tute for the pressure but it gives an acceptable fit to the simulations 
discussed below. 

Given M, a/M , M, ap, T, Y, a nd r B , the slim disk equations 
provide q(r), and me lGammiel d 19991) equations provide b- r (r), b^{r), 
and p(r). 

The remaining free parameters are a t and n. These parameters 
can be inferred from MHD simulations. Note that oq and n control 
the size and shape of the turbulent contribution to a(r), and ai and 
r B control the size and shape of the mean field contribution to a(r). 

In § 94161 we estimate values for these four parameters by 
fitting the a(r) prescription to data from six GRMHD simula- 
tions. As an example, Figure [TJ shows a{r), as prescribed by equa- 
tion l |30t , for the parameters of Model A in Table [4] (solid red 
curve). The mean magnetic field contribution to a(r) (long-dashed 
green curve) dominates inside the ISCO, where the plunging fluid 
stretches and amplifies the magnetic field. The turbulent contribu- 
tion to a(r) (dashed blue curve) dominates outside the ISCO, where 
mean magnetic fields are weak. At large radii, the shear parame- 
ter becomes 3/2 and a becomes constant. Data from a GRMHD 
simulation (gray points; c.f. 95. 51 are in good agreement with the 
one-dimensional prescription for a(r). 

In the next three sections, we detail our GRMHD simulations 
and their connection to the a(r) prescription. 



4 DETAILS OF THE SIMULATIONS 
4.1 Computational Method 

The simulations were carried o ut wi th the 3D G RMHD 
code HARM dGammie et al.l 120031 ; iMcKinnevI 120061 : 



iMcKinnev & Blandfordl |2009|) . which solves the ideal MHD 
equations for the motion of a magnetized gas in the Kerr metric, 
the spacetime of a rotating black hole. The equation of motion 
of the gas is taken to be u = p/(T - 1), where u and p are the 
internal energy and pressure and T is the adiabatic index. The code 
conserves energy to machine precision, so any energy lost at the 
grid scale by, e.g., turbulent dissipation or numerical reconnection, 
is returned to the gas, increasing its entropy. 

TableQ]gives a summary of the six simulations, which we have 
labeled A-F. Simulations A, B, and C are thin, radiatively efficient 
disks, and simulations D, E, and F are thick, radiatively inefficient 
disks. The spin parameter is a/M = 0.7 for simulations B and E, 
and a/M = for the others. 

Most of the simulations have been described in previous pa- 
pers, so our overview of the simulations in this section can be 
brief. Simulations B and C are two of the models discussed by 
Kulkar ni et al. I d201lh . Simulations A, B, and C were analyzed by 



Zhu et al. (where they are labeled E, B , and C, respec t ively) . 



Finally, simulations D and F were studied bv lNaravan et al.l |2012|) 
(where they are called SANE and MAD). The only simulation that 
has not appeared before is E, but it only differs from D in that it has 
a/M = 0.7 and a duration of 100, 000M. 

The resolution of simulations B and C in r x 9 x <f> is 256 x 64 x 
32, and the resolution of the other simulations is 256 x 128 x 64. 
The radial grid is logarithmically spaced to concentrate attention 
on the inner regions of the flow. The inner boundary of the grid is 
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Figure 1. The a(r) prescription denned by equation j30t . for parameters 
oq = 0.025, a\ = 100, n = 6, and rg = 6M (solid red curve). This pre- 
scription is the sum of two terms, a mean magnetic field component (long- 
dashed green), which dominates inside the ISCO, and a turbulent compo- 
nent (dashed blue), which dominates outside the ISC O. These two com- 
ponents are base d on the one-dimensional models of Gammlel 1 19991) and 
Sadowski (201 1), respectively. At large radii, a(r) converges to ag = 0.025, 
a constant. This model is a good description of the data from simulation 
A (gray points; c.f. ^5}. Data from inside r str j cl are marked with filled cir- 
cles, data from between r slr j ct and r\ oost , are marked with open circles, and 



data from outside ri oosc are marked with crosses (c.f. 35.5 



between the Cauchy horizon and event horizon, and outflow bound- 
ary conditions are used there, so the event horizon behaves as a true 
horizon. The polar grid is squeezed towards the equatorial plane to 
concentrate resolution on the turbulent, high density regions of the 
flow, at the expense of the laminar, coronal re gions. Simulations D, 
E, and F use a version of the grid developed bv lTchekhovskov et al.l 
d201 in which the 6 resolution near the pole increases with in- 
creasing radius so as to follow the formation of jets, which colli- 
mate at large distances. The azimuthal grid is uniform and extends 
from to max , where max is either tt/2 (simulations A, B, and C) 
or 2n (simulations D, E, and F). 

The properties of the six simulations are summarized in Table 



1 



4.2 Initial Conditions 

Initially, the gas orbits the black hole in a torus in hydro static 
equilibrium dDe Villiers et alj2003l : |Penna et alj201(ll2012al) . The 
thickness of the torus can be adjusted to give either thin or thick ac- 
cretion disks. A weak poloidal magnetic field threads the torus. All 
of the simulations start with a sequence of poloidal loops, except F, 
which starts from a single magnetic loop. When there are multiple 
loops, the black hole accretes flux of alternating polarity over time 
and little net flux builds up on the hole. In simulation F, the center 
of the loop at r = 300M does not reach the black hole over the du- 
ration of the simulation, so the black hole acquires a large net flux. 
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Table 1. GRMHD simulation parameters 



Simulation 


a/M 


h/r 


Initial loops 


n r xn g x n$ 


0max 


Duration 


A 





0.1 


Multiple 


256 x 128 x 64 


71 


20, 000M 


B 


0.7 


0.05 


Multiple 


256 x 64 x 32 


nil 


27, OOOM 


C 





0.05 


Multiple 


256 x 64 x 32 


Jc/2 


27, OOOM 


D 





0.3 


Multiple 


256 x 128 x 64 


In 


200, OOOM 


E 


0.7 


0.3 


Multiple 


256 x 128 x 64 


In 


100, OOOM 


F 





0.3 


Single 


264 x 126 x 60 


In 


100, OOOM 



In all of the simulations, the magnetic field is normalized so that 
the initial gas-to-magnetic pressure ratio has minimum /? = 100. 

4.3 Quasi-Steady State 

The initial condition is unstable to the MRI. Differential rotation of 
the torus triggers the MRI and the gas becomes turbulent after ~ 10 
orbits. Turbulence transports angular momentum and energy out- 
wards and the gas accretes inwards. At late times in the simulation, 
gas in the disk near the midplane is turbulent. Magnetic buoyancy 
lifts fields above and below the disk, forming a highly magnetized 
corona. The corona is mostly laminar because the MRI requires 
p>\. 

Figure [2] shows the fluid frame magnetic field at the end of 
simulations A and D. The field has been azimuthally averaged but 
not time averaged. The coordinates are x/M = rsin(#) and z/M = 
rcos(0). The turbulent region of simulation A is thinner than the 
turbulent region of simulation D. The turbulent region of simulation 
D extends nearly to the polar axes. 



5 ANALYSIS OF SIMULATIONS A AND D 

In this section we discuss our analysis of simulations A and D. 
These describe a prototypical thin disk and a prototypical thick disk 
around non-spinning black holes. We discuss the remaining four 
simulations in §6\ 

Our goal for this section, achieved in 35.61 is to extract a{r) 
profiles from the simulation data and compare them with the one- 
dimensional prescription defined by equation l |30| >. As a first step, 
we discuss the distinction between disk and coronal fluid. We only 
include disk fluid in our calculations. Then we discuss the radial 
range of fluid that can be considered to have reached a quasi-steady 
state. We only include quasi-steady data in our calculations. We 
examine the shear rate and epicyclic frequency of the simulations, 
because these play an important role in determining a. Finally, we 
compute a(r) and compare it with our prescription from ^3] 

5.1 The Distinction Between Disk and Corona 

We would like to separate the disk component of the flow from 
the coronal component so that we can focus our analysis on the 
disk. There are several reasons to do this. For one, the stress has 
a different character in the corona and disk regions of the flow. In 
the corona, the stress is mostly generated by mean magnetic fields, 
while in the disk, the stress is generated largely by turbulence. So 
including coronal stresses in the model would add new difficulties. 
For this reason, and also for simplicity, we focus on the disk region 
of the flow. 



There are other reasons to isolate the disk from the corona. At 
least in thin accretion disks, the emission from the corona and disk 
are different. The disk has a thermal spectrum and the corona has a 
power law spectrum. So the distinction is sensible for observations. 
Accretion disk models which use the a- viscosity prescription tend 
to focus on the disk region of the flow and ignore the corona. An- 
other reason to separate out the disk region is that our numerical 
grid concentrates 6 resolution at the midplane and leaves the polar 
regions poorly resolved. So the simulation data are unreliable in the 
corona. 

We therefore only include fluid within one density scale height 
of the midplane in our analysis. The density scale height is defined 
as, 

where p is rest mass density in the fluid frame, and pu' is rest 
mass density in Boyer-Lindquist coordinates. The time integral 
is over the steady state period of the flow, as explained below. 
Another popular definition for the scale height is (h/r) rms = 
(j (O-n/2) 2 p^gdtd9d</}/ f p^[=gdtded<p) ' '. We have no rea- 
son to favor one definition over another, though it sho uld be noted 
that ( fe/y)rm« can be a factor of ~ 2 bigger than h/r l lPenna et al.l 
l201Ch . 

We also add a density-weighting to vertical averages, to fur- 
ther emphasize midplane fluid. That is, the density weighted verti- 
cal average of O is jOpu' -\fgmd9/ f pu' sfg^dO. 

The top panel of Figure[3]shows log(a) as a function of x and 
z for simulation A. We explain how a is obtained from the simula- 
tion in 35.51 but we show this plot here to illustrate the distinction 
between the disk and the corona. Black dashed curves mark one 
scale height above the midplane. Note that a has a very different 
character above and below the disk region. In the coronal regions 
a is much larger than in the disk. The bottom panel of Figure [3] 
shows the ratio of Maxwell stress to Reynolds stress on a loga- 
rithmic scale. Maxwell stresses are much more significant in the 
corona. The different character of the flow in these two regions is 
one of the reasons we only include disk fluid within one h/r of the 
midplane in our calculations. 

Figure|4]shows the same quantities for simulation D. Again, a 
and the Maxwell stress are much larger outside the disk. However, 
the shape of this region does not track h/r as it did in simulation 
A. In fact, the high a, high Maxwell stress region has a parabolic 
shape, bounded by roughly z/M = (r/6M) . This looks like a jet. 
For simplicity and consistency, we will restrict our calculations to 
the fluid within one h/r of the midplane for all of the simulations. 
We have checked that our results do not depend on the details of 
this cutoff, as long as we do not include the "jet" region. 
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averaged to make the streamlines appear continuous in this two-dimensional projection (this makes the flow appear slightly less turbulent). Turbulent twisting 
of the magnetic field can be seen on different scales. When the field is twisted on the smallest scale, the grid scale, there is reconnection and dissipation. The 
density scale height of the disk is h/r ~ 0. 1 (c.f. 35. U . Fluid in the disk region of the flow is turbulent. Magnetic buoyancy lifts magnetic fields out of the disk 
where they settle in a highly magnetized coronal region. The coronal region is mostly laminar, because magnetic tension is quenching the MRI. On very large 
scales, the magnetic field has an approximately dipolar structure. Right panel: Same as left panel, but for simulation D. This accretion flow is much thicker 
and the flow only becomes laminar near the polar axes. 



5.2 Radial and Azimuthal Velocities 

To extract smooth results from turbulent data, it is necessary to av- 
erage the data over at least a viscous tim e, which smooths out tur- 
bulent fluctuations dNaravan et alj|2012h . As a first step, we need 
the radial velocity of the simulations. To compute the shear rate 
and epicyclic frequency of the gas (and thus a), we need the az- 
imuthal velocity of the gas. We compute these two components of 
the velocity in this section. We will not need the 6 component of 
the velocity. It is much smaller than the radial and azimuthal com- 
ponents in the disk region of the flows. 

In a thin accretion disk, the gas closely follows circular 
geodesies as it spirals toward the ISCO. The radial velocity of 
the gas is significantly smaller than the azimuthal velocity. In a 
standard thin disk, the radial velocity is suppressed by a factor of 
a(h/r) 2 . In the ADAF solution, which describes a very thick accre- 
tion flow, the radial velocity is suppressed by a factor of a relative 
to the azimuthal velocity. These relations hold approximately for 
the GRMHD simulations as well. 

The radial velocity most relevant for turbulence is the one 
measured by the ze ro angular momentum observer (ZAMO) of 
iBardeen et al.l dl972h . This is a local, inertial frame attached to ob- 
servers with zero angular momentum. As a result of frame drag- 
ging, ZAMO observers appear to rotate with respect to observers at 
infinity. The ZAMO frame is at rest with respect to the local space- 
time. 

The radial velocity in the ZAMO frame is 



VaV 
~7F 



(32) 



function of radius for simulation A. The data has been averaged 
over the disk region of the flow, with a density weighting, as dis- 
cussed in 95.11 It has been time-averaged over the quasi-steady 
state period of the flow, which we explain below. The ISCO is at 
r = 6M. Inside the ISCO, the gas is approximately in free fall and 
the radial velocity increases rapidly as the gas approaches the black 
hole. Outside the ISCO, the motion is more nearly circular. In this 
region, the radial velocity is suppressed relative to the azimuthal 
velocity as predicted by standard disk theory. Because the radial 
velocity is small, it takes a long time for the simulation to reach a 
quasi-steady state. Once steady state is reached, it requires a time 
average extending over many orbital periods to smooth out turbu- 
lent fluctuations and obtain reliable results. We take up these issues 
in the next section. 

The top right panel of Figure [5] shows the radial velocity as 
a function of radius for simulation D. The radial velocity is much 
larger than the radial velocity of A, because this accretion flow is 
geometrically thick. For this reason, the data from this simulation 
is in quasi-steady state out to a larger radius. Also, the larger radial 
velocity smears out the inner edge of the accretion disk. There is 
no longer a sudden transition between slow and fast radial velocity, 
at the ISCO or at any other radius. 

The bottom left panel of Figure[5]shows the angular velocity 



«= — = —, 

dt u' 



(33) 



where A and A are defined by equations ([2j and (|4j. 

The top left panel of Figure [5] shows the radial velocity as a 



of simulation A as a function of radius. Again, we have in- 
cluded only fluid within one scale height of the midplane, taken 
the density-weighted vertical average, and time averaged over the 
steady state portion of the flow. The dashed line shows the an- 
gular velocity of circular geodesies in the equatorial plane, D = 
M i/2 /(r il2 + aM 1 ' 2 ). The flow follows circular geodesies except for 
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x/M 



Figure 3. Top panel: log(ti) in the r — 8 plane for simulation A. The data has 
been time-averaged over t = 7, 000M - 20, 000M. Dashed black lines indi- 
cate one density scale height above and below the midplane (c.f. 35. i\ . We 
refer to the low a region within one scale height of the midplane as the disk, 
and the high a region outside one scale height as the corona. We restrict our 
calculations to the disk, for the reasons discussed in 35.11 We show this 
plot here to illustrate the difference between the disk and the corona. We 
explain how a is obtained from the simulation in 35.51 Bottom panel: Time- 
averaged ratio of Maxwell stress to Reynolds stress on a log scale, in the 
r - 8 plane, for simulation A. Coronal fluid is more magnetically dominated 
than disk fluid. This is expected, as magnetic buoyancy lifts magnetic fields 
into the corona. 

r < 5M, where the radial velocity is increasing rapidly. The angu- 
lar velocity peaks around 3M. At this radius, the shear parameter q 
must go to zero, so the turbulent contribution to a becomes negli- 
gible. 

The bottom right panel of Figure [5] shows the angular veloc- 
ity of simulation D. It is also nearly Keplerian outside the ISCO. 
This is perhaps surprising, because the angular velocity of the self- 
similar ADAF solution is very sub-Keplerian. The initial torus of 
the GRMHD simulation persists at large radii over the duration 
of the simulation and continues to feed nearly Keplerian gas into 
the inner regions of the flow. This acts as a very strong boundary 
condition, which may limit the solution's ability to converge to the 
ADAF solution. The ADAF solution is self-similar and describes an 
accretion flow with infinite extent but, as we will see, the simulation 
is only converged out to r = 100M. A longer duration simulation, 
which has reached quasi-steady state out to a larger radius, might be 




50 100 150 200 

x/M 




50 100 150 200 



x/M 

Figure 4. Same as Figure [5] except for simulation D. Unlike simulation A, 
the highly magnetized region does not have the same shape as the disk scale 
height. In fact, the former has a paraboloidal, jet-like shape, approximately 
z/M = (r/6M) 2 . 

expected to have a more ADAF-like angular velocity. Nonetheless, 
the radial velocity of simulation D does appear to be co nverging to 
the ADAF prediction, as shown bv lNaravan et al.l d2012b . 

5.3 Convergence and Steady State 

Following iNaravan et alj d2012l) , we divide the data from simula- 
tion D into six "time chunks" which are logarithmically spaced 
in time. Each time chunk is about twice as long as the previous 
one. They are summarized in Table [3] This logarithmic spacing 
is useful since most of the quantities we are interested in show 
power-law behavior as a function of both time and radius. Note 
that there is no overlap between chunks, and hence each chunk 
provides independent information. Because the duration of simu- 
lation A is only 20, 000M, we use a single time chunk, spanning 
t = 7, 000M - 20, 000M. 

For each time chunk, we compute the time-averaged ra- 
dial velocity profile v,(r) of the gas within one scale-height 
of the mid-plane. We estimate the viscous time at radius r by 
dNovikov &Thornell 1973k IPenna et aI1l2010h : 

fvispW = \ — TTT' (34) 

\v r (r)\ 

Following INaravan et al.l ( |2012|) . we then define two criteria, one 



© 0000 RAS, MNRAS 000, 000-000 



10 R. F. Penna, A. S. Sqdowski, A. K. Kulkarni, & R. Narayan 




Figure 5. Top left: Radial velocity as a function of radius for simulation A. The solid line extends to r = r str j ct and the dashed line extends to r = noose 
(estimated convergence radii of 35.31 . The radial velocity increases suddenly around the ISCO at r = 6M, inside of which there are no stable circular orbits 
for the gas to follow. Top right: Radial velocity as a function of radius for simulation D. Colors correspond to time chunks 1 (blue), 2 (green), 3 (red), 4 
(cyan), 5 (magenta), and 6 (black) (see 35.31 . Bottom left: Angular velocity as a function of radius for simulation A. The dashed black curve shows the angular 
velocity of circular, equatorial geodesies. The simulated flow is nearly geodesic outside the ISCO. The angular velocity has a maximum near the photon orbit 
at r = 3M, so the shear parameter (and hence the turbulent contribution to a) will be zero here. Bottom right: Angular velocity as a function of radius for 
simulation D. This flow is slightly sub-Keplerian outside the ISCO. 



"strict" and one "loose," to estimate the radius range over which 
the flow has achieved inflow equilibrium: 

'vise (^strict) = 'chunk/2, (35) 
'vise (/loose) = 'chunk, (36) 

where f cnU nk is the duration of the chunk. The values of f cnun k, ''strict, 
and r loose for the various time chunks are summarized in Tables [2] 
and [3] We only trust data from inside ri oose . Data from inside r str j Ct 
are considered particularly reliable. 



5.4 Shear and Epicyclic Frequencies 

Now we compute the dimensionless shear parameter, q, and the 
epicyclic frequency, k, as a function of radius. The shear parameter 
is computed from q = -Txr^jQ.. The epicyclic frequency is com- 
puted from the shear parameter using equation ( 125b . The data are 
vertically averaged over the gas within one scale height of the mid- 
plane using the density weighting. The data are time averaged over 
the time chunks listed in Tables|2]and[3] 

The top left panel of Figure [6] shows the shear parameter as a 
function of radius for simulation A out to r strict (solid curve) and 
loose (dotted curve). The analytical shear parameter for circular, 
equatorial Kerr geodesies is shown for comparison (dashed curve). 
At large radii, the analytical q converges to the shear parameter 
of non-relativistic Keplerian flow, q = 3/2. At the ISCO, general 
relativistic corrections increase the shear parameter to q = 2. The 
GRMHD shear parameter is about 10% larger than the analytical 
shear parameter outside the ISCO. Inside the ISCO, the analytical 
q blows up as it approaches the photon orbit. The GRMHD shear 



parameter goes to zero near the photon orbit and is negative very 
close to the black hole. 

The top right panel of Figure[6]shows the shear parameter as a 
function of radius for simulation D. Results are shown for each of 
the time chunks. All of the time chunks are consistent out to ri oose to 
within several percent. This gives us confidence that the simulation 
has converged to a quasi-steady solution. The GRMHD shear pa- 
rameter is similar to the shear parameter of simulation A. It is about 
10% larger than the analytical q outside the ISCO, turns over inside 
the ISCO, and drops to zero near the photon orbit. There is good 
agreement between the GRMHD and analytical shear parameters 
out to r loose ~ 100M. 

The bottom left panel of Figure [6] shows the epicyclic fre- 
quency as a function of radius for simulation A. The bottom right 
panel shows the same quantity for simulation D. In both cases, out- 
side the ISCO, the epicyclic frequency of the simulation is about 
10% lower than the epicyclic frequency of Keplerian flow. In both 
cases the epicyclic frequency has a minimum near the ISCO. The 
minimum of the epicyclic frequency roughly marks the most unsta- 
ble radius in the flow, because k = corresponds to marginal stabil- 
ity. So it is consistent with standard disk theory that the minimum 
of k is near the ISCO. Simulation D has a broader and shallower 
minimum, indicating the inner edge of this disk has been "smeared 
out" by the larger radial velocity of the flow. 



5.5 Shakura-Sunyaev Viscosity Parameter, a 

Finally, we compute the dimensionless viscosity parameter, a. The 
GRMHD stress-energy tensor is a combination of Reynolds and 
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Table 2. Convergence radii for simulations A, B, and C 



Simulation 


Time Range (M) 


'chunk /M 


r s uict/M 


loose/M 


A 


7,000-20,000 


13,000 


9 


10 


B 


20,000-27,000 


7,000 


6.5 


7 


C 


20,000-27,000 


7,000 


6.5 


7 



Table 3. Time chunks for simulation D 



Chunk 


Time Range (M) 


'chunk IM 


r s tnct/M 


Ooose IM 


I 


3,000-6,000 


3,000 


19 


23 


II 


6,000-12,000 


6,000 


25 


43 


III 


12,000-25,000 


13,000 


29 


45 


IV 


25,000-50,000 


25,000 


43 


62 


V 


50,000-100,000 


50,000 


66 


92 


VI 


100,000-200,000 


100,000 


86 


113 



Maxwell terms: 

T, y = + T<r g) , (37) 

where, 

T^v y) = (P + + pV> ( 38 ) 

r<r E) = ^ (&Vv + bl K* - . ( 39 ) 

and = + u^iiy is the projection tensor. To each term there 
is an associated stress, which is the rip component of the tensor 
measured in the fluid frame. So the Reynolds stress is and the 

Maxwell stress is T^ ag> . An important difference between the two 
is that the Reynolds stress requires turbulence whereas the Maxwell 
stress can be generated by turbulence or large scale magnetic fields 
and thus can be nonzero even in laminar flow. 

We define a as the ratio of total stress to total pressure: 



We have chosen to include b 2 /2 in the denominator because it keeps 
a < 1 . As a result of this choice, part of the observed dependence of 
a on q is inherited via the magnetic pressure because the magnetic 
pressure is amplified by shear. This is significant inside the ISCO 
where the magnetic pressure is comparable to or exceeds the gas 
pressure. 

We compute the Maxwell and Reynolds stresses in the rest 
frame of the mean flow, 

1 f*"" 

S" = if dip, (41) 

<pmax JO 

rather than in the rest frame of the instantaneous flow, it. The mean 
flow is the 0-averaged instantaneous flow. There are some subtleties 
in the distinction between these two frames. The Reynolds stress 
vanishes in the rest frame of the instantaneous flow, u p , by defini- 
tion of the fluid frame, but not in the rest frame of the mean flow, 
uf. This is our primary reason for choosing uf to define the iner- 
tial fluid frame when computing a . The electric field vanishes in 
the rest frame of the instantaneous flow, by the assumption of ideal 
MHD, and not in the rest frame of the mean flow, but for simplicity 
we assume the electric field can be neglected in both frames. 



One could inc lude a time av erage in the definition of the mean 
flow, equation {41}. Penna et al. (2010) included a 100M time aver- 
age. That is, they averaged the instantaneous flow over all of <p and 
over 100M in time to obtain the mean flow. At large radii, 100M 
is much smaller than the orbital timescale, so this extra averaging 
has little effect. However, inside the ISCO, 100M is larger than the 
orbital timescale. In this case, the time-averaged mean flow tends 
to give larger a. The reason is that time-averaging increases the 
discrepancy between the mean and instantaneous flows by adding 
contributions to the mean flow from earlier and later times. This 
discrepancy propagates into a when the stress tensor is boosted to 
the rest frame of the mean flow. 

Figure|7]shows a(r) for simulation A with and without a 100M 
time averaging in the definition of W . Including the time averaging 
increases a and the effect is greatest inside the ISCO. In fact, the 
peak a exceeds unity when the mean flow is defined with a time 
average. This would imply the stresses carry more energy than the 
total energy of the gas and magnetic fields, which is unphysical. 
To avoid these sorts of contradictions, we do not include any time 
averaging in equation <4U . 

The middle panel of Figure [8] shows the ratio of Maxwell 
stress to Reynolds stress in simulation A. Outside the ISCO, 
the ratio tracks the prediction of the linearized MRI, (4 - q)/q 
dPessah et alj|2006bh . In fact, it is slightly larger, as is consistent 
with s hearing box simulations of MRI turbulence jPessah et al.l 
2006b). Inside the ISCO, the Maxwell stress is an order of mag- 
nitude larger than the Reynolds stress. The flow is mostly lami- 
nar inside the ISCO, so Reynolds stress is weak, and the plunging 
fluid stretches the magnetic field, so the Maxwell stress is strong. 
The bottom panel of Figure [8] shows the pr oduct a/j. It approaches 
the expected value of » 0.5 in the disk (Blackman et al J [20081 : 
iGuan et alj2009l : ISorathia et alj|2012l) . 

The top panel of Figure [9] shows a(r) for the six time chunks 
of simulation D. Outside the ISCO, the various time chunks agree 
to within 30%, which provides an estimate for the contribution of 
turbulent noise to the error in our a measurements. Inside the ISCO, 
where the flow is more nearly laminar, the agreement between the 
time chunks is better. 

The bottom panel of Figure [9] shows the ratio of Maxwell to 
Reynolds stress as a function of radius for simulation D. Inside 
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Figure 6. Top left: Dimensionless shear parameter, q, as a function of radius for simulation A (solid and dotted curves) and for Keplerian flow (dashed). The 
simulation data are plotted out to r\ oox (dotted curves) and out to r str j ct (solid curves). The simulated shear parameter turns over near the ISCO. Top right: 
Dimensionless shear parameter as a function of radius for simulation D. Colors are as in Figure|5] Bottom left: Epicyclic frequency as a function of radius for 
simulation A. The epicyclic frequency has its minimum near the ISCO, as expected. Bottom right: Epicyclic frequency as a function of radius for simulation 
D. The minimum is still near the ISCO, but it is broader and shallower than the minimum in the epicyclic frequency of simulation A. The larger radial velocity 
of simulation D has "smeared out" the inner edge of the disk. 



about r = 22M, the ratio is consistent with the ratio of simulation 
A. Outside r = 22M, the ratio begins to grow with radius. It is not 
clear what causes this. It may indicate that simulation D has not 
reached quasi-steady state at these radii yet. Figure [4] shows highly 
magnetized, irregular clumps of fluid in the disk region, even after 
time-averaging the data over the last time chunk. In a true quasi- 
steady state, one would expect time averaging to eliminate these 
clumps. 



5.6 Comparison With the One-Dimensional a(r) Prescription 

Finally, we compare a{r) of the simulations with the one- 
dimensional prescription for a(r) of ^3] as defined by equation i30\ 
and the parameters listed in Table|4] 

Figure[T]shows the agreement between the GRMHD a(r) from 
simulation A and the a(r) prescription. Figure [10] shows the agree- 
ment between simulation D and the a{r) prescription. The inter- 
pretation of the GRMHD a(r) in terms of a mean magnetic field 
component in the inner regions, and a turbulent component in the 
outer regions, appears to match the data. It is more difficult to fit 
simulation D because the a{r) prescription relies on a sharp distinc- 
tion between the laminar, magnetically dominated inner regions of 
the flow, and the turbulent, weakly magnetized outer regions. This 
distinction is cleanest in thin disks, like simulation A, which have 
a clear inner edge at the ISCO. In thick disks, like simulation D, 
the inner edge of the disk is smeared out by the radial velocity of 
the flow, so the separation of a into two independent components 
is less sound. 

The shape of a{r) does not depend strongly on all of the pa- 
rameters in Table [4] The mass of the black hole in the GRMHD 
simulations is dimensionless. It is listed in Table |4]because it goes 



into the slim disk solutions that underlie the one-dimensional a{r) 
prescription. We set M/M a = 10 arbitrarily and this choice has a 
negligible effect on a(r). 

The accretion rate also enters through the slim disk part of 
a(r). Our estimates of M/M e dd for th e GRMHD simulations are 
based on the analysis of IZhu et"al1 ( l2012h . who used h/r as a proxy 
for the accretion rate. This gives rough estimates, which is all we 
need because the dependence of a(r) on the accretion rate is also 
very weak. 

The magnetic flux threading the horizon, T, is measured di- 
rectly from the GRMHD simulations. It slightly affects the shape 
of a(r) inside the ISCO. 

The four parameters that strongly control the shape of a(r) are 
a , a,, n and r B . It is encouraging that a = 0.025, and n = 6 give 
good fits to both simulations. In other words, both simulations have 
a(r) ~ 0.025 at large radii, and a cc q 6 (and q > 0) in the turbulent 
disk. 

Mean magnetic fields are only important inside the ISCO of 
simulation A, so we set r B = 6M in this case. The region where 
mean magnetic fields are important in simulation D is broader, so 
we set r B = 30M in this case. 



6 ANALYSIS OF SIMULATIONS B, C, E, AND F 

In this section we present data from four more GRMHD simula- 
tions. This gives information about the dependence of the viscosity 
parameters a?o, ati, n, and r B , on black hole spin, resolution, and 
the amount of magnetic flux threading the black hole. Of these ef- 
fects, the dependence on flux threading the black hole is the most 
dramatic. 
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Table 4. Parameters for a(r) fits to the GRMHD simulations 



Simulation 


M/Mq 


a/M 


M/Medd 


T 


ao 


Of] 


/i 


i~b 


A 


10 





0.5 


0.6 


0.025 


100 


6 


nsco 


B 


10 


0.7 


0.2 


3 


0.025 


10 


6 


Isco 


C 


10 





0.2 


6 


0.025 


1 


6 


Isco 


D 


10 





1 


5 


0.025 


0.5 


6 


30M 


E 


10 


0.7 


1 


10 


0.025 


0.5 


6 


30M 


F 


10 





1 


30 


0.025 


0.1 


6 


30M 




0.01 

— 1000 

It 100 

> 10 

E* > 

0.1 
0.01 



Simulation A 



H 1 1 1 1 1 1 — |- 



-H ' — I 1 1 1 1 — I — |- 




Figure 7. Dimensionless viscosity parameter, a, as a function of radius for 
simulation A, with (magenta) and without (black) a 100M time-average 
in the definition of the mean fluid frame, h' 1 . Outside the ISCO, 100M is 
longer than an orbital period, so the effect is small. Inside the ISCO, 100M 
is several orbital periods, so the effect is large. 



6.1 Simulation E 

Simulation E is identical to simulation D except the black hole has 
spin parameter a/M = 0.7 and the duration is 100, 000M. We have 
divided the simulation data into time chunks as we did for simu- 
lation D, but there is one less time chunk because the duration is 
half as long. The time chunks and radii of convergence estimates 
are listed in Table [5] 

Figure fTTIshows the radial and angular velocities as a function 
of radius for the five time chunks. After the first three time chunks, 
around t = 25,000, the radial velocity drops by about a factor of 
2. Then it holds steady (to within a few percent) for the final two 
time chunks. So even at t = 25, 000M, the simulation is still settling 
down. The radial velocity at the end of simulation E is about 30% 
lower than the radial velocity at the end of simulation D. As a result, 
simulation E has converged over a smaller range of radii; we find 
loose = 60M for the last time chunk (simulation D had r\ oox = 90M 
over the same time interval.) 

Figure [T2] shows the dimensionless shear parameter and 
epicyclic frequency as a function of radius. Outside the ISCO, the 
shear parameter is about 20% larger than the shear parameter of cir- 



Figure 8. Top panel: Dimensionless viscosity parameter, a, as a func- 
tion of radius for simulation A. The data has been time averaged from 
t = 7, 000M to 20, 000M. Solid and dotted curves correspond to v ^ ^strict 
and r < Hoosei Middle panel: Time-averaged ratio of Maxwell to Reynolds 
stress for simulation A (solid and dotted curves). The dashe d curve is the 
the pre diction from a linearized MRI analysis, A/q(r) — 1 IPessah et alj 
2006b), for Keplerian q(r), equation 1241 . Bottom panel: Product aB for 
simulation A (solid and dotted curves). The dash ed line at aB = 0.5 is 
the expected value for saturated MR I turbulence (Blackm an et al J [20081 : 
iGuan et al.l l2009; Sorathia et al. 2012). The simulated product falls below 
this value inside the ISCO where the flow is mostly laminar. 



cular geodesies and the epicyclic frequency is about 20% smaller, 
similar to the results from simulation D. The minimum epicyclic 
frequencies of the two simulations are also comparable. Simulation 
D had k/Q. ~ 0.6 and simulation E has k/Q. ~ 0.5. 

Figure[T3]compares a(r) as computed from the last time chunk 
of the simulation against a(r) computed from the one-dimensional 
viscosity prescription with a = 0.025, q-j = 0.5, n = 6, and r B = 
30M. This is the same choice of parameters that gave a good fit 
to simulation D. It is interesting that they fit simulation E as well. 
This suggests the parameters of the modified viscosity prescription 
do not depend strongly on black hole spin. 
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Table 5. Time chunks for simulation E 



Chunk 


Time Range (M) 


fchunk/-^ 


''strict IM 


noose IM 


I 


3,000-6,000 


3,000 


9.5 


15 


II 


6,000-12,000 


6,000 


15 


20 


III 


12,000-25,000 


13,000 


22 


31 


IV 


25,000-50,000 


25,000 


31 


44 


V 


50,000-100,000 


50,000 


44 


60 
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6.2 Simulation B 

This simulation is identical to simulation A, except the black hole 
is spinning with spin parameter a/M = 0.7 and the resolution is 
256 x 64 x 32 rather than 256 x 128 x 64. 

The left panel of Figure [JJ] shows the GRMHD a(r). The a(r) 
prescription is shown for the same parameters that gave a good fit 
to simulation A: a = 0.025, ati = 1, n = 6, and r B = r lsC o- The 
parameters oq = 0.025 and n = 6 governing the turbulent part of 
a(r) are the same across all four simulations we have considered 
so far, suggesting these parameters do not depend strongly on disk 
thickness. 

6.3 Simulation C 

Simulation C is the same as simulation A except the resolution is 
lower: 256 x 64 x 32 versus 256 x 128 x 64 and the disk is thinner 
(h/r ~ 0.05 instead of h/r ~ 0.1). The data from this simulation 
has a lower a than the a(r) prescription with our fiducial choice 
of parameters a = 0.025 and n = 6. This suggests we are under- 
resolving the MRI at this resolution. To infer whether the values 
ao = 0.025 and n = 6 are robust, higher resolution simulations 
would be useful. 

The results are shown in the right panel of Figure [74l 

6.4 Simulation F 

Simulation F differs from simulation D in one crucial respect. The 
initial torus of gas is threaded with a single poloidal magnetic field 
loop rather than multiple loops. The center of the initial loop is 
centered at r = 300M and gas from this radius does not reach the 
black hole over the duration of the simulation. So the polarity of 
the flux that reaches the black hole is approximately const ant and 
a large net flux builds up on the hole. Nar avan et alj J2012h give a 
detailed account of the convergence in time and radius, and the role 
of outflows, in simulations D and F. 

The large net flux carried by the gas in simulation D has a 



dramatic effect: the flow remains mostly laminar at all radii. Figure 
[731 shows the fluid frame magnetic field in the r — plane at t = 
100, 000M, the final snapshot of the simulation. The eddies and 
turbulent twisting of the field are all but gone on every scale, in 
marked contrast to the other five simulations we considered (see, 
e.g., Figure|2j. 

Following lNaravan et al.|j2012h . we divide the simulation data 
into five time chunks. The time periods and estimated convergence 
radii for each time chunk are summarized in Table [6] This simu- 
lation has the largest radial velocity of any of the simulations (see 
Figure [T6"ll. so the estimated convergence radii are the largest. The 
final time chunk has r strict = 170M and r ]oax = 207M. 

Simulation F also has the most sub-Keplerian angular velocity 
of the six simulations (Figure [TBI. The angular velocity drops by a 
factor of a few between time chunks I and III, but it is consistent 
across the final three time chunks to within a few percent. 

The ratio of Maxwell stress to Reynolds stress has a much 
clumpier distribution in the r — 8 plane than any of the other simu- 
lations. A large, magnetized, z-shaped clump, where the Maxwell 
stress is enhanced, extends throughout the flow (bottom panel of 
Figure [TTt. The irregular shape of the clump suggests it is a non- 
equilibrium structure. Perhaps if the duration of the simulation was 
longer it would be smoothed out. The appearance of the magne- 
tized clump suggests r strict = 170M is a better estimate for the radii 
of convergence for this simulation than ri oose = 207. 

The shear parameter of the flow (top panel of Figure|18l> varies 
by about 50% between time chunks. Outside the ISCO, the shear 
parameter is roughly constant with radius. The epicyclic frequency 
(bottom panel of Figure [T8t has its minimum near r = 20M, rather 
than at the ISCO. The minimum itself is very broad and shallow, 
not extending much below k/Q. = 1. In other words, the inner edge 
of the disk has moved well outside the ISCO and is highly smeared 
out. 

The profiles of a as a function of radius for the five time 
chunks are shown in the top panel of Figure [19] Outside r « 20M, 
the profiles of a are constant with radius, even increasing slightly. 
The other simulations have a decreasing with radius. This suggests 
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Figure 14. Same as Figure[T]except for simulations B (left panel) and C (right panel). 



Table 6. Time chunks for simulation F 



Chunk 


Time Range (M) 


fchunk/M 


raneilM 


loose /M 


I 


3,000-6,000 


3,000 


35 


52 


II 


6,000-12,000 


6,000 


37 


65 


III 


12,000-25,000 


13,000 


69 


90 


IV 


25,000-50,000 


25,000 


109 


128 


V 


50,000-100,000 


50,000 


170 


207 



the turbulent contribution to a is not dominating even at the largest 
converged radii, which is consistent with the laminar structure of 
the magnetic field lines. 

The bottom panel of Figure [T9] shows the ratio of Maxwell to 
Reynolds stress as a function of radius for the five time chunks. 
For the first two time chunks, the Maxwell stress is an order of 
magnitude larger than the Reynolds stress at all radii. During the 
final three time chunks, the ratio appears to have stabilized out to 
r = 40M. At larger radii, the Maxwell stress is enhanced by the 
non-equilibrium, magnetized, z-shaped clump noted earlier. 



Our a(r) prescription does not appear to give a good fit to 
the simulation F results (Figure l20t. This is probably because the 
simulation is mostly laminar at all radii (as shown by Figure [T5ll, 
whereas our a(r) prescription assumes turbulence dominates the 
stress beyond the innermost radii. For simulations A-E, this is a 
good assumption provided the disk region of the flow is distin- 
guished from the coronal regions. However, the entire domain of 
simulation F is mostly laminar and highly magnetized, and so it 
should perhaps be considered entirely coronal gas. It appears a dif- 
ferent a(r) prescription is needed to describe such highly magne- 
tized flows. 



7 a-DISK SOLUTIONS WITH VARIABLE a(r) 

In this section, we evaluate the dependence of a-disk solutions 
on the a(r) prescription. The particular q-disks we c onsider are 
"slim disks" dAbramowicz et alj[l988l ; ISadowskillioTH) . We com- 
pare slim disk solutions with constant a = a Q to solutions with 
varying a(r), where a(r) is defined by equation d30t and the param- 
eters inferred from simulation A (c.f. Tabled- That is, we consider 
a non-spinning, 10 solar mass black hole, threaded with a magnetic 
flux T = 0.6. The viscosity parameters are ao = 0.025, n = 6, 
ai\ = 100, and r B = 6M. We consider two different accretion rates: 
30% Eddington and Eddington. 

Figure [2T1 shows our results. At large radii, a(r) converges to 
a , so the solutions with constant and varying a(r) are the same 
to within a percent. Inside the ISCO, the fluid plunges toward the 
black hole with little dissipation, so in the innermost regions the 
solutions are again insensitive to the a prescription. Only in an in- 
termediate zone, between the ISCO and r « 20M, does the shape 
of a(r) have a significant effect. For solar mass black holes, this 
region emits predo minately in X-ray s and is relevant for black hole 
spin measurements dGou et al ,11201 lh , 

In this zone, the a(r) prescription has a larger a than the con- 
stant a = a prescription. So the a{r) prescription increases the 
disk's radial velocity by a factor of 2 - 3 and lowers its central 
(mid-plane) temperature by about 50%. In fact, a solution accret- 
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Figure 16. Same as Figure \E\ but for simulation F. Of the six simulations 
we consider, this simulation has the largest radial velocity and is the most 
sub-Keplerian. 



ing at the Eddington limit with varying a{r) has the same central 
temperature as a solution accreting at 30% Eddington with constant 
a (c.f. Figure[2Tt. So the a(r) prescription has a significant effect on 
central temperature. Central temperature depends on both effective 
temperature and optical depth, so the effect is really due to changes 
in surface density. 

At low accretion rates, the disk is radiatively efficient and the 
effect of a(r) on the radiated fluxed is negligible. At high accretion 
rates, advection becomes important and the radiated flux shows its 
dependence on a(r). At the Eddington limit, the flux from the so- 
lution with varying a(r) is about 50% lower than the flux from the 
solution with constant a. So flux is only affected by the a(r) pre- 
scription at high accretion rates. 



8 DISCUSSION AND SUMMARY 

The a(r) prescription of ^3] must be computed numerically. How- 
ever, for all practical purposes, the function a(r) for thin disks can 
be reduced to a simple analytical formula. Our simulations suggest 



a(r) = 0.025 



q(r) 



3/2 



(<? > 0), 



(42) 



where q{r) is given analytically by equation I l24t . The constant co- 
efficient and exponent in equation {42} are the values favored by 
our GRMHD simulations (c.f. Table They will change as bet- 
ter simulation data becomes available. The free parameters are the 
mass and spin of the black hole, which enter through equation J24b 
for q{r). The contribution from mean magnetic fields can be ig- 
nored in this approximation because mean field stresses are only 
significant inside the ISCO, and thin disks are not sensitive to a(r) 
in this region. Equations ( 124b and J42b thus give an analytical a(r) 
prescription that can be used for thin disk models (see Figure [22l>. 
The more general a(r) prescription of ^3]is needed for thick disks. 
To summarize our main results, we have constructed a 





™o 



100 

x/M 



Figure 17. Same as Figure [5] but for simulation F. A large, highly- 
magnetized, z-shaped clump persists over most of the flow inside ri oose . 
The magnetized clump increases a and the ratio of Maxwell to Reynolds 
stress. 



one-dimensional prescription for a(r) and estimated parameters 
for this prescription based on data from GRMHD simulations. 
The f act that a varies w i th radius had been antic ipated long 
ago jPringle & Reeslli972l ; Ishakura & Sunvaevlll973h but global 
MHD simulations provide the first quantitative information about 
the sh ape of a(r) jPapaloizou & Nelsonl 120031 : IPenna et alj|20ld 
l2012bMFromang et al]201lh . 

Our modified a{r) prescription, equation j30K is the sum of 
two components. The first component describes mean field stresses. 
It is important in the laminar, inner regions of accretion disks, 
where the plunging fluid stretches the frozen-in field. Our descrip - 
tion of this component is based on the model of iGammiej l i 1999b . 
which supplies the magnetic stress as a function of radius and two 
free parameters: black hole spin and the amount of magnetic flux 
threading the horizon. 

The second component of th e a(r) prescription de scribes tur- 
bulent stresses. As emphasized bv lPessah et ai] ]2006b), a depends 
on the shear parameter, q. In Newtonian gravity, Keplerian flow has 
a constant shear parameter, q = 3/2, but general relativistic correc- 
tions give even Keplerian disks around black holes a varying q(r), 
as discussed in Sj2] The shear parameter increases from q = 3/2 at 
the outer edge of the disk to q = 2 at the inner edge of the disk. 
This is a 50% change in q but it creates a larger variation in a, be- 
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Figure 19. Same as Figure[8]but for simulation F. 



cause a goes as q" (for q > 0). Our GRMHD simulations are too 
noisy to infer n precisely, but the data seem to prefer n x 6. This 
is consistent with the simulations of Pess ah et al. (2006b), which 
resulted in n between 2 and 8. Analytical MHD closure models for 
the MRI also allow n between 2 and 8 (Kato & Yoshizawj[l99l 
19951: lQgilviell2003l : iPessah et alj|2006j|bll2008h . 

Simple extensions of the standard a pre scription and some 
closure models predict ne gative a for q < jKato & Yoshizaw3 
1993lll995l:[Ogil vie 2003). An exception is the closure model of 
Pessah et al . (2006a b). Data from shearing box simulations show 
zero turbulent stress for q < JPessah e t al. 2008). Our simulations 




10 15 zo 



r/M 



Figure 21. Slim disk solutions with varying a(r) (solid curves) and with 
constant a = 0.025 (dashed curves), for the parameters inferred from simu- 
lation A (c.f.[4). Solutions are shown at the Eddington accretion rate (blue 
curves) and at 30% Eddington (red curves). Top panel: Radial velocity as 
a function of radius. Solutions with varying a(r) have larger radial veloc- 
ities. Middle panel: Midplane temperature (not effective temperature) as a 
function radius. Solutions with varying a(r) are colder in the X-ray emit- 
ting portions of the flow. Bottom panel: Radiant flux as a function of radius. 
Radiant flux is affected by a(r) only at high accretion rates. 
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Figure 22. Analytical a(r) prescription defined by equation J421 for a/M = 
(solid red curve). This model gives a good fit to data from Simulation A. 
Data from inside r stl -; ct are marked with filled circles, data from between 
r strict an d loose s[e marked with open circles, and data from outside ri 00se 
are marked with crosses (c.f. 35. 51 . 



are consistent with this data, although they are not as decisive on 
this point because negative shear parameters are only found inside 
roughly the photon orbit, where mean field stresses are large. Our 
prescription for a(r) is always positive and the turbulent contribu- 
tion vanishes for q < 0. The mean field term in our a(r) prescription 
gives a good description of our simulation data in regions near the 
black hole where q < 0. 

We have discussed accretion onto black holes. When a disk 
accretes onto a star, a boundary layer forms between the star and 
disk. It can generate half the accretion luminosity in soft X-rays 
JPringld[l977h . The boundary layer in stellar accretion is similar 
to the region inside the ISCO in black hole accretion. In both re- 
gions, the angular veloc ity is non-Keplerian and the shear ampli- 
fies t h e magnetic field (Armitage 2002; St einacker & Papaloizou 
l2002h . lSteinacker & Papaloizod feoolv found a(r) profiles in MHD 
simulations of boundary layers that are similar to our a(r) profiles 
inside the ISCO. It would be interesting to extend the a(r) prescrip- 
tion to these cases. 

We analyzed data from six GRMHD accretion disk simula- 
tions. Three of the simulations are thin, radiatively efficient accre- 
tion disks (simulations A, B, and C). The other three are geomet- 
rically thick, radiatively inefficient accretion flows (simulations D, 
E, and F). The simulations vary in resolution from 256 x 64 x 32 
to 256 x 128 x 64. Two of the simulations describe spinning black 
holes, with spin parameter a/M = 0.7 (simulations B and E) and 
the others describe non-spinning black holes. MRI driven turbu- 
lence and large scale magnetic fields generate stresses in the sim- 
ulated disks self-consistently, so the a viscosity prescription is not 
assumed. Instead, we measure a(r) from the simulation data. 

For each simulation, we measure a{r) = T- r i/p tot and compare 
it against the a(r) prescription. Thin accretion disks are easier to de- 
scribe with the a(r) prescription than thick accretion disks because 
they have a sharp transition at the ISCO that creates two distinct 



regions: a magnetically dominated region inside the ISCO, and a 
weakly magnetized, turbulent region outside the ISCO. This dis- 
tinction is blurred in thick accretion disks by the large radial veloc- 
ity of the flow. So some of the simplifications in the a(r) prescrip- 
tion are less applicable. Simulation F, in which a large magnetic 
flux was allowed to build up on the hole, is particularly difficult to 
interpret because turbulence is almost completely absent. 

We were careful to only include fluid which has reached a 
quasi-steady state. The timescale to reach quasi-steady state scales 
as the inverse of radial velocity, and is thus an increasing function 
of radius. So the inner regions of disks converge before the outer re- 
gions. Thin accretion disks have smaller radial velocities than thick 
accretion disks, so thin disk simulations are only converged out to 
r ~ 10A/ rather than r ss 100M. 

For further insight into the simulations, we analyzed their 
shear parameters and epicyclic frequencies. Outside the ISCO, the 
shear parameters of the simulations are usually within 20% of the 
Keplerian prediction. Inside the ISCO, the shear parameter turns 
over, going to zero near the photon orbit at r = 3M where the 
angular velocity peaks. So the turbulent contribution to a, which 
scales as q" (for q > 0), is unimportant near the black hole. 

The epicyclic frequency of the flow is also close to the Kep- 
lerian value outside the ISCO. The inner edge of the flow can be 
identified with the minimum in K(r), the radius where the flow is 
most unstable. This is usually near the ISCO, although when the 
disk is thick the inner edge is smeared out by the large radial veloc- 
ity of the flow. 

Finally, we considered the effect of the a(r) prescription on 
cy-disk solutions. We compared solutions with varying a(r) to solu- 
tions with constant a = Oq. We fixed the free parameters ao,ai,n, 
and r B using the values inferred from simulation A (c.f. Table [4]l. 
The differences between varying and constant a are only significant 
in the region between the ISCO and r as 20M. At smaller radii, the 
gas is plunging too quickly for stresses to act, so a does not en- 
ter, and at larger radii the varying a(r) prescription is converging to 
a = cto- In the intermediate zone between the ISCO and r « 20M, 
the solutions with varying a(r) have larger a than the solutions with 
constant a. This increases their radial velocity and lowers their cen- 
tral temperature and radiant flux. The effect on central temperature 
is most significant. A solution accreting at the Eddington rate with 
radially varying a(r) has the same central temperature as a solution 
accreting at 30% Eddington with constant a. Central (mid-plane) 
temperature depends on both the effective temperature and the op- 
tical depth, so the effect is really due to changes in surface density. 
The effect of a on the radiant flux is only important at high accre- 
tion rates. 

The main shortcomings of our a(r) prescription are the ab- 
sence of gas pressure in the prescription for the mean magnetic field 
component, and the absence of magnetic fields in the prescription 
for the turbulent component. For thin, weakly magnetized disks 
these are better assumptions than for thick or highly magnetized 
disks because the two components do not overlap significantly in 
the disk. 

The simulations are limited by their duration, which prevents a 
large range of radii from reaching quasi-steady state. They are also 
limited by their resolution. The fastest growing mode of the MRI is 
usual l y not resolved by mo re than about 10 grid cells jPenna et al.l 
l20ld:]Naravan et alj|2012h which is only marginally acceptable 



Shiokawa et alj|2012l ; ISorathia et alj|2012h . Despite these limita- 



tions, these are among the best GRMHD accretion disk simulations 
available at present. Shearing box simulations can resolve the local 
physics of the MRI better, but cannot obtain the dependence of a 
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on radius explicitl y. Nonetheless, the shearing box simulations of 
IPessah et alj J2008I) should be revisited with higher resolution (they 
used 32 x 192 x 32), as these are the best way to measure the de- 
pendence of a on q. 

Our results suggest that relativistic corrections to q partly 
contribute to the higher q 's measured in GRMHD simulations 
( [Penna et alj l20ldL l2012bl) compared to Newto nian simulations 
(Papaloizou & Nelson! 2003: Fro mang et alj|201lh . Assuming a cc 
q 6 (for q > 0), we infer that a is six times larger in the relativis- 
tic inner regions of GRMHD disks than in Newtonian disks. This 
partly resolves the discrepancy but does not go far enough, as the 
a's measured from GRMHD simulations are over an order of mag- 
nitude larger than the a's measured from Newtonian MHD simula- 
tions. More work is needed to understand this discrepancy. 

Switching from a constant a to a radially v arying a(r) would 
have a small effect on black hole spin estimates. iGou et al] ilOl lh 
considered the effect on continuum fitting measurements of the spin 
of the black hole in Cyg X-l if one assumes a = 0.01 versus 
a = 0. 1 . The black hole spin decreases slightly, from a/M = 0.9988 
to 0.9985, as a is increased. Switching from a constant a to the 
a(r) prescription will have a similar effect. This is well below the 
current observational sources of error in bla ck hole spin measure- 
ments dKulkarni et al.ll201lUZhu et alj2 012). so it is not a concern 
for now. 

Black hole spin estimates are restricted to observations of 
disks radiating below 30% of the Eddington limit , which corre- 
sponds to thin disks (h/r ~ O.l. lKulkarni et alj201lh . Observations 
based on models of thick disks will be more sensitive to the shape 
of a(f). 
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